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PREFACE 


The  work  presented  in  this  technical  report  was  a  joint  effort 
between  researchers  at  the  Air  Force  Research  Laboratory  ( AFRL) 
and  the  University  of  Florida.  The  Computational  Mechanics 
Branch  (AFRL/MNAC)  was  the  primary  developer  of  the  high  order 
ENO  solution  algorithm  for  the  governing  equations  of  uniaxial 
strain  (ID)  and  plane  strain  (2D) .  The  University  of  Florida 
developed  and  implemented  the  interface  resolution  and  tracking 
algorithms . 


1.  Introduction 


The  dynamics  of  high-speed  impact  between  materials  is  characterized  by  large 
deformation  and  short  time  scales.  Wave  propagation  in  the  impacting  media  is 
highly  nonlinear,  and  involves  localized  phenomena  such  as  shear  bands,  crack 
propagation,  and  wave  refraction  (Meyers,  1994).  These  problems  are  typically 
challenging  to  solve  because,  in  contrast  to  conventional  structural  dynamics 
problems,  the  deviatoric  and  pressure  terms  in  the  stress  tensors  are  both  important 
and  need  to  be  modeled  separately.  In  contrast  to  conventional  fluid  dynamics 
problems,  the  stress  and  strain  fields  are  related  through  nonlinear  elasto-plastic  yield 
surfaces,  the  models  for  which  must  be  included  in  the  governing  equations. 
Furthermore,  the  interface  between  materials  experiences  not  only  fast  motion,  but 
also  large  variations  in  shape.  Such  problems  have  been  routinely  handled  by  the  so- 
called  hydrocodes.  Benson  (1992)  provides  a  very  detailed  review  of  the  formulation, 
modeling  and  computational  techniques  employed  by  these  large-scale  computer 
codes. 

In  the  past  several  years,  progress  has  been  made  toward  developing  predictive 
techniques  based  on  the  application  of  modem  computational  solid  and  fluid 
dynamics.  For  example,  Camacho  and  Ortiz  (1996,  1997)  have  developed  a 
Lagrangian  finite  element  impact  dynamics  model  involving  brittle  materials 
(Camacho  and  Ortiz  1996),  and  ductile  penetration  (Camacho  and  Ortiz  1997).  Their 
approach  is  based  on  adaptive  meshing,  explicit  contact/friction  algorithm,  and  rate 
dependent  plasticity.  Trangenstein  (1990,  1994,  1995) ,  and  Trangenstein  and  Pember 
(1991)  have  adopted  Godunov’s  method  and  ideas  developed  in  modern 
computational  fluid  dynamics  to  handle  the  problem  in  the  context  of  a  Riemann  type 
problem  with  second-order  accuracy.  In  addition,  there  have  also  been  finite  volume 
methods  based  on  general  discretization  treatments  (Bailey  and  Cross  1995,  Dormy 
and  Tarantola  1995)  with  first-order  accuracy. 

Overall,  it  seems  clear  that  while  progress  has  been  made,  impact  dynamics  has  been 
rigorously  investigated  only  recently.  Much  remains  to  be  done  before  we  can  treat  it 
satisfactorily  in  terms  of  geometric  representation,  interface  movement,  temporal  and 
spatial  resolution,  accuracy  of  numerical  schemes  and  incorporation  of  realistic 
materials  behavior  models. 

In  this  report  we  describe  the  development  of  a  numerical  solution  technique  for  the 
simulation  of  high-speed  multimaterial  impact.  Of  particular  interest  is  the  interaction 
of  solid  impactors  with  targets.  This  problem  is  important  in  applications  such  as 
munition-target  interactions,  geological  impact  dynamics,  shock  processing,  and 
formation  of  shaped  charges  upon  detonation  and  their  subsequent  interaction  with 
targets  (Meyers,  1994).  Such  interactions  present  the  following  challenges  to 
numerical  simulation  techniques: 
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1 .  High  velocities  of  impact  leading  to  large  deformations  of  the  impactor  as  well  as 
targets. 

2.  Nonlinear  wave-propagation  and  the  development  of  shocks  in  the  systems. 

3.  Modeling  of  the  constitutive  properties  of  materials  under  intense  impact 
conditions  and  accurate  numerical  calculation  of  the  elastic -plastic  behavior 
described  by  the  models. 

4.  Phenomena  at  multiple  interfaces  (such  as  impactor-target,  target-ambient  and 
impactor-ambient),  i.e.  both  free  surface  and  surface-surface  dynamics. 

This  report  describes  a  numerical  approach  that  seeks  to  tackle  all  of  the  above 
aspects.  In  the  following  description  we  elaborate  on  the  numerical  framework  for 
solution  of  high-velocity  impact  problems  and  present  a  rationale  for  choice  of  the 
numerical  methodology.  We  also  address  issues  of  accuracy  in  the  treatment  of  the 
wave  propagation  in  the  bulk  media  and  interactions  at  the  interfaces,  and  present 
results  showing  the  capability  of  the  method. 

2.  Governing  Equations  for  ID  case:  Uniaxial  strain 

The  system  being  simulated  in  the  ID  case  is  as  shown  in  Figure  1.  We  consider  two 
slender  rods  impacting  on  the  x-axis.  The  rods  can  be  of  finite  or  infinite  extent.  The 
system  of  governing  equations  for  the  elastic-plastic  flow  under  uniaxial  strain 
conditions  is  of  the  form: 
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Figure  1.  Impact  of  two  bars  along  the  x-axis.  Setup  for  ID  uniaxial  strain  case. 


where  p  is  the  density,  u  is  the  velocity,  E  is  energy,  e  is  the  equivalent  plastic  strain, 
sx  is  the  deviatoric  stress. 

The  flux  vector  is  of  the  form 
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are  the  source  terms  arising  in  the  equations  for  deviatoric  stress  and  equivalent 
plastic  strain.  G  is  the  shear  modulus  of  the  material  and  ae  =  H(e  )  is  the  equivalent 

stress-strain  curve  as  determined  by  uniaxial  strain  experiments.  The  pressure  is 
determined  by  an  equation  of  state  (eos)  of  the  form: 


p  =  eos(p,pu2,E ) 


(7) 


The  specific  eos  in  use  here  is  the  Mie-Gruneisen  eos  (Meyers  1994): 
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where  po  is  the  density  of  the  unstressed  material,  Co  and  s  are  coefficients  that  relate 
the  shock  speed  Us  and  the  particle  velocity: 


Us  =c0  +  su 


(12) 


The  conditions  that  apply  at  the  ends  of  the  rods  depend  on  the  extent  of  the  rods  and 
are  addressed  later  in  this  report. 


3.  Governing  equations  for  the  2D  case:  Plane  strain 

In  2D  the  impact  of  arbitrary  shapes  can  be  considered  in  the  plane  as  illustrated  in 
Figure  2.  The  equations  in  2D  correspond  to  the  plane-strain  situation  and  are  as 
follows: 


The  transport  equation  in  vector  form  is: 


dQ  {  dF(Q)  |  dG(Q)  =  d 
dt  dx  dy 


(13) 
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Figure  2.  Impact  of  two  objects  in  the  2D  case.  A  plane  strain  problem  is  solved. 


The  vector  of  independent  variables,  and  the  x-  and  y-direction  convective  flux 
vectors  are  given  by: 


(14b) 
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The  source  vector  is: 
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where  the  source  term  in  the  equation  for  equivalent  plastic  strain  is: 
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In  this  and  other  expressions  where  the  tensor  index  notation  is  mixed  with  Cartesian 
component  notation,  subscripts  1,2  or  3  correspond  to  x,y  or  z  respectively.  The 
source  term  in  the  equation  for  deviatoric  stress  component  is: 
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In  the  above  the  rate  of  deformation  tensor  components  D,j  are  used  and  these  are 
defined  as 
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where  8,j  is  the  Kronecker  delta.  The  terms  in  Equation  (16b)  that  involve  Q  are 
needed  to  maintain  material  frame  indifference. 

In  the  equation  set  above  the  three  components  of  deviatoric  stress  are  chosen  as  the 
dependent  variables.  The  total  stress  components  are  related  to  the  deviatoric  stresses 
as: 

°ij=sij-p8ij  (18) 


where  by  definition  p  =  — — okk . 

In  addition,  one  again  requires  the  eos  for  the  pressure  which  is  given  by  the  Mie- 
Gruneisen  equation  of  state,  Eq.  (8).  The  above  set  of  equations  is  solved  in  the 
presence  of  moving  boundaries,  namely  the  impactor  and  target.  The  numerical 
approach  adopted  for  solving  the  governing  equations  needs  to  incorporate  the 
presence  of  these  moving  boundaries  where  various  types  of  boundary  conditions 
may  apply.  Therefore  the  critical  issue  of  the  manner  in  which  interfacial  dynamics  is 
approached  numerically  is  addressed  next. 

4.  Numerical  Approaches  for  Moving  Boundary  Problems  in  Fluid  Mechanics 

A  large  variety  of  physical  phenomena  involve  the  coupling  of  evolution  of  flowfields 
or  material  deformation  fields  with  boundaries  that  move,  deform  or  evolve  in  time. 
Examples  include  the  deformation  of  drops,  bubbles,  liquid  free  surfaces  etc.,  the 
evolution  of  phase  boundaries  in  solidification  and  vaporization,  fluid-structure 
interaction  problems  at  the  large  scale  such  as  in  aeroelasticity,  and  in  the  small  scale 
such  as  in  biomechanics  and  a  whole  host  of  other  interesting  phenomena.  These 
problems  are  challenging  to  CFD  practitioners  due  to  the  complexity  associated  with 
the  often  severely  deformed  boundaries  in  (or  enclosing)  the  flow,  and  the 
nonlinearity  resulting  from  the  coupling  of  the  interface  dynamics  with  the  dynamics 
of  the  flowfield.  Ideally  one  would  like  to  track  the  moving  boundary  as  a  sharp  front 
(assuming  discontinuities  in  some  flow  quantity  such  as  stress,  energy  etc.  are  present 
across  the  interface)  without  smearing  information  at  the  front.  Also,  one  would  like 
to  solve  the  flowfield  within  the  regions  separated  by  the  interfaces  with  desired 
accuracy.  If  the  interfaces  become  multiply  connected  one  would  wish  to  follow  the 
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evolution  of  the  interfaces  through  such  topological  changes.  Of  course,  it  would  also 
be  desirable  to  correctly  model  the  physics  so  that  all  the  ingredients  of  the  dynamics 
can  be  included  in  the  model.  Algorithms  for  computational  solution  of  moving 
boundary  problems  operate  under  these  demands. 

In  the  computational  fluid  dynamics  literature  numerous  methods  have  been 
developed  for  the  purpose  of  handling  moving  boundary  interaction  with  flows.  These 
can  be  classed  broadly  under  the  categories  of  Eulerian,  Lagrangian  and  mixed 
(Eulerian-Lagranian)  methods.  We  will  briefly  summarize  the  salient  features  of  each 
method  in  the  following: 

1.  Lagrangian  methods:  Ideally  one  would  want  to  simulate  the  effect  of  the 
boundary  by  treating  it  explicitly,  without  smearing  the  information  at  the 
interface,  i.e.  with  minimum  numerical  diffusion.  There  are  several  ways  of  doing 
this.  For  a  fixed  boundary,  when  the  shape  is  truly  complex,  one  can  resort  to 
block-structured  domain  decomposition  (Shyy,  1994),  overset  meshes  (Steger, 
1999,  Johnson  and  Belk,  1995),  or  unstructured  boundary-conforming  curvilinear 
grids  (Venkatakrishnan,  1996)  to  discretize  the  domain.  The  last  method  has  been 
developed  very  strongly  by  the  finite-element  community  (Camacho  and  Ortiz, 
1996,  1997)  and  has  been  used  also  by  the  fluid  dynamics  community  (Fritz  and 
Boris,  1979,  Glimm  et  al.,  1988,  Rausch  et  al.,  1993)  for  moving  boundary 
problems.  There  has  also  been  an  emergence  of  meshless  methods  for  the  solution 
of  finite  deformation  problems  (Liszka  et  al.,  1996,  Duarte  et  al.,  1996,  Randles  et 
al.,  1996).  These  have  been  successful  in  solving  an  impressive  array  of  problems. 
For  moving  boundaries  which  may  undergo  large  deformations,  or  are  subject  to 
topology  changes  in  the  course  of  their  evolution,  generating  body-fitted  grids  to 
conform  to  complex  boundaries  may  become  taxing.  Furthermore,  as  will  be 
elaborated  below  there  are  certain  advantages  to  solving  large  distortion  problems 
on  fixed  structured  meshes  in  terms  of  straightforward  extensions  of  well- 
developed  numerical  schemes,  application  of  fast  solvers,  absence  of  grid 
distortion,  entanglement  and  massive  re-gridding  effects,  and  the  ability  to  follow 
the  evolution  of  interfaces  through  changes  in  topology  without  affecting  the 
computational  grid.  There  are  of  course  some  significant  issues  involved  in  making 
fixed  grid  methods  work  and  these  are  also  mentioned  below. 

2.  Eulerian  methods:  In  this  class  of  methods,  particularly  for  surfaces  undergoing 
large  deformations  such  as  in  free-surface  flows  or  shattering  of  droplets,  it  is 
advantageous  to  dispense  with  tracking  the  interface  as  a  curve  or  surface. 
Instead,  in  this  class  of  methods  the  boundary  is  deduced  from  values  on  the  mesh 
of  a  field  variable  which  could  be  a  volume-of-fluid  (Hirt  and  Nichols,  1981, 
Brackbill  et  al.,  1992,  Kothe  and  Mjolsness,  1992,  Scardovelli  and  Zaleski,  1999), 
a  level-set  (Osher  and  Sethian,  1988),  or  phase-field  (Kobayashi,  1993,  Wheeler 
et  al.,  1992),  or  the  enthalpy  in  solidification  problems  (Voller  and  Prakash, 
1987).  The  interface  is  then  an  isocontour  of  the  appropriate  field  variable. 
Usually  calculations  are  performed  over  a  fixed  Cartesian  mesh.  However, 
recently  some  Eulerian  methods  have  been  extended  to  curvilinear  (Zhang  et  al.. 
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1998),  adaptive  (Sussman,  1997)  and  unstructured  meshes  (Barth  and  Sethian, 

1998) .  Eulerian  methods  perform  very  well  for  a  variety  of  moving  boundary 
problems  and  are  perhaps  ideally  suited  to  free-surface  problems.  However,  in 
these  problems,  particularly  when  surface  forces  are  to  be  included  in  the  flow 
calculations,  the  interface  essentially  is  diffuse  and  occupies  a  few  grid  cells  in 
practical  calculations.  This  may  be  undesirable  in  many  problems  both  from  an 
accuracy  as  well  as  physical  reality/modeling  standpoint.  The  problem  treated  in 
the  current  report  is  one  example. 

3.  Mixed  (Eulerian-Lagrangian)  methods:  This  class  of  methods  combines  some 
advantages  of  the  two  approaches  above.  The  interface  is  tracked  explicitly  as 
curves  (or  surfaces  in  3D).  The  computations  are  performed  on  fixed  meshes 
whose  topology  is  independent  of  that  of  the  interface.  One  example  of  this  type 
of  method  is  the  immersed  boundary  technique  used  for  a  range  of  multifluid 
problems  (Unverdi  and  Tryggvason,  1992,  Juric  and  Tryggvason,  1996, 
Udaykumar  et  al.,  1997)  and  for  fluid-structure  interaction  problems  in  biofluids 
(Peskin,  1977,  Fauci  and  Peskin,  1988).  While  explicitly  tracking  the  interface, 
this  method  transmits  the  information  regarding  the  discontinuity  across  the 
interface  to  the  grid  in  much  the  same  way  as  purely  Eulerian  methods,  i.e.  by 
casting  the  surface  forces  into  a  body  force  term  in  the  governing  equations. 
Therefore  the  solution  reverts  to  a  one-domain  approach,  i.e.  the  solver  does  not 
see  a  discontinuity  at  the  location  of  the  interface,  but  experiences  distributed 
forces  and  material  properties  in  the  vicinity  of  the  interface.  As  demonstrated  by 
Beyer  and  Leveque  (1992)  this  results  in  a  method  that  is  globally  0(h)  accurate, 
where  h  is  the  grid  spacing.  On  the  other  hand,  the  cut-cell  treatment  (Udaykumar 
and  Shyy,  1995b,  Udaykumar  et  al.,  1996,  Udaykumar  et  al.,  1999,  Ye  et  al., 

1999)  proceeds  to  reconstruct  the  domain  on  either  side  of  the  interface  with 
attention  to  the  immersed  boundary  and  its  geometry  overlying  the  grid.  Phases 
are  treated  separately  and  no  smearing  of  the  interface  takes  place  at  the 
formulation  level.  A  conservative  control  volume  treatment  demanding  care  in 
assembly  of  fluxes  is  performed  for  cells  on  both  sides  of  the  interface.  This 
method  has  been  applied  in  recent  work  to  problems  involving  immersed 
stationary  and  moving  solid  boundaries  in  incompressible  flows  where  global 
second-order  accuracy  has  been  demonstrated. 

The  choice  of  moving  boundary  method  from  the  general  categories  above  depends  to 
a  large  extent  on  its  appropriateness  for  the  physical  problem  chosen.  In  this  report 
we  demonstrate  a  fixed  grid  approach  for  high-speed  impact  that  can  be  a  powerful 
tool  for  numerical  computation.  The  method  presented  in  this  report  falls  under  the 
class  of  mixed  methods.  It  operates  on  a  fixed  Cartesian  mesh  (the  Eulerian  part) 
while  the  interfaces  move  through  the  mesh  (the  Lagrangian  part).  The  method  treats 
the  interfaces  as  discontinuities  without  smearing  on  the  mesh,  therefore  it  is  a  sharp 
interface  method.  The  advantage  of  the  fixed  grid  approach  is  obviously  that  grid 
topology  remains  simple  while  large  distortions  of  the  interface  take  place.  The  issues 
involved  in  explicit  interface  tracking  are  discussed  below. 
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When  interfaces  are  tracked  explicitly,  periodic  reorganization  of  the  interface 
information  becomes  necessary.  This  can  result  from  dilation  or  compression  of 
segments  of  the  interface  or  due  to  topology  changes  of  the  interface.  In  2D,  mergers 
and  breakups  can  be  handled  quite  effectively  (Udaykumar  and  Shyy  1995a,  Juric  and 
Tryggvason  1996).  In  3D  the  operations  can  become  more  complicated  (Snyder  and 
Woodbury,  1993).  Therefore,  in  3D  situations,  explicit  tracking  of  interfaces  will  be 
work-intensive  in  the  context  of  mergers  and  breakups  in  comparison  to  purely 
Eulerian  methods.  However,  there  are  physical  problems  for  which  explicit  interface 
information  becomes  desirable;  then  interface  tracking  becomes  the  natural  choice  for 
immersed  boundary  treatment.  One  instance  is  when  a  solid-liquid  boundary  is  being 
tracked  and  where  the  no-slip  boundary  condition  is  to  be  applied.  This  can  be  done 
for  boundary-fitted  grid  computations  quite  naturally.  For  fixed  grid  computations 
also  this  can  be  done  if  the  exact  location  of  the  interface,  as  provided  by  explicit 
tracking,  is  known  (Udaykumar  and  Shyy  1997,  Udaykumar  et  al.  1999,  Ye  et  al. 
1999).  In  fluid-structure  interaction  problems,  such  as  in  the  dynamics  of  membranes 
(Dong  et  al.,  1988,  Kan  et  al.,  1998),  adhesion  of  cells  to  substrates  in  biofluids 
(Jones  et  al.,  1995)  or  the  dynamics  of  pliable  aerodynamic  surfaces  (Shyy  et  al. 
1996;  Smith  and  Shyy,  1996;  Fauci  and  Peskin,  1988),  the  forces  generated  within  the 
membranes  depend  on  the  stretching  and  bending  of  the  membranes.  This  requires 
information  on  the  dilatation  of  the  interface.  Furthermore,  for  the  case  where  the 
membrane  is  anchored  to  a  surface  as  in  the  adhesion  of  a  cell  membrane  to  a 
substrate,  the  forces  transmitted  to  the  membrane  need  to  be  calculated  (Dembo  et  al., 
1988).  Explicit  tracking  is  ideally  suited  to  these  types  of  situations.  The  ability  of 
mixed  Eulerian-Lagrangian  methods  to  incorporate  both  solid-liquid  no-slip 
boundaries  as  well  as  fluid-fluid  interfaces  has  been  demonstrated  in  previous  work 
by  the  authors  (Udaykumar  et  al.,  1997). 

In  the  problem  tackled  in  this  report,  an  Eulerian-Lagrangian  method  can  be  used  to 
good  effect.  Here  we  wish  to  simulate  the  interaction  of  a  solid  impactor  with  targets 
and  to  track  the  interface  separating  the  two  materials  in  time  to  determine  the 
deformation  of  the  two  interacting  bodies.  When  the  materials  are  in  contact 
following  impact  the  material  discontinuities  at  the  interface  need  to  be  tracked 
without  numerical  diffusion.  Explicit  interface  tracking  is  ideally  suited  to  this  task. 
Having  chosen  to  explicitly  track  the  boundaries,  a  fixed  grid  is  used  to  compute  the 
flow  solutions  in  the  presence  of  the  moving  boundaries.  This  allows  an  extension  of 
highly  accurate  shock-capturing  methods  such  as  Essentially  Non-Oscillatory  or  ENO 
(Harten  et  al.,  1997,  1987,  Shu  and  Osher,  1988,  1989)  in  the  present  case.  These 
methods  have  been  developed  for  scalar  conservation  laws  in  fixed  grid  settings  to 
solve  moving  boundary  problems  with  arbitrarily  distorted  interfaces.  A  Cartesian 
grid  ENO  formulation  suffers  little  change  when  applied  to  the  present  problem.  We 
now  proceed  to  describe  the  method  in  detail. 
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The  interface  is  described  by  interfacial  markers  defined  by  the  coordinates  Xfsj.The 
spacing  between  the  markers  is  maintained  at  some  fraction  of  the  grid  spacing, 
0.5h<ds<1.5h.  The  convention  adopted  is  that  as  one  traverses  the  interface  along  the 
arclength,  the  material  enclosed  by  the  interface  lies  to  the  right.  This  is  illustrated  in 

Figure  3.  The  functions  x(s)  =ax  s?  +  bx  s  +  cx  and  y(s)  =a-y  +  by  s  +  c-y  are 
generated.  The  coefficients  <^x/y,>  b^y  and  cx/y  at  any  interfacial  point  i  are  obtained 
by  fitting  polynomials  through  the  coordinates  (*/./, y;.;),  ( x/,y/, )  and  (xi+j,yi+j). 

The  coefficients  ax/y>t  b^y  and  cx/y  are  stored  for  each  marker  point.  The  normal  to 
the  interface  then  points  from  the  interior  to  the  exterior  of  the  object  and  is  given  by 
the  equation: 


n  =  ( 


-yx 


yR  +  f,2  Vxv2  +  ^-2 


(19) 


The  derivatives  xs,  ys  are  evaluated  using  central-differencing  along  the  arclength 

coordinate  s.  Cubic  splines  were  also  tried  without  observable  differences  in  the 
results  for  previous  test  problems.  Therefore,  central  differencing  was  adopted  since  it 
is  easily  applicable  to  different  end  conditions  for  the  boundaries. 
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Figure  4.  Interfacial  cell  and  bulk  cell  classification  on  a  grid  with  an  interface 
passing  through  it.  Also  shown  are  interfacial  cell  properties. 


6.  Relationship  between  interface  and  grid 

Once  the  interface  has  been  defined,  the  information  on  its  relationship  with  the  grid 
has  to  be  established.  There  may  be  several  interfaces  (henceforth  called  objects) 
immersed  in  the  domain.  Each  of  the  objects  may  enclose  material  with  different 
transport  properties.  Therefore  it  is  necessary  to  identify  which  phase  each 
computational  point  (i.e.  cell  center  point)  lies  in.  The  procedure  for  obtaining  this 
and  related  information  has  been  discussed  in  detail  in  Udaykumar  et  al.(1999).  The 
end  result  of  the  procedures  are  the  following  pieces  of  information  which  are 
required  to  set  up  the  discretization  scheme  for  the  present  method  : 

1 .  The  interfacial  cell  in  which  each  interface  marker  lies. 

2.  The  interfacial  marker  which  is  closest  in  distance  to  a  computational  point. 

3.  The  material  in  which  each  computational  point  in  the  mesh  lies. 

4.  Several  geometric  details  such  as  the  shape  of  the  resulting  cut-cell,  the 
locations  where  the  interface  cuts  the  cell  faces  and  where  it  intersects  the  cell 
center  lines  (the  dotted  lines  shown  in  Figure  4).  These  details  of  a  cell  are 
used  in  constructing  the  stencil  for  each  interfacial  cell. 

5.  A  list  of  all  interfacial  cells. 

These  pieces  of  information  regarding  the  interface  and  its  relationship  to  the 
underlying  grid  are  computed  only  in  a  lower-dimensional  set  of  interface  cells. 
Therefore  using  local  searches  and  operations  and  data  storage  limited  to  this  set  of 
cells  renders  dealing  with  the  interface  and  mesh  relationship  economical.  In  practical 
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runs  the  operations  associated  with  obtaining  the  interface  and  mesh  information  is  a 
small  fraction  of  the  computing  time  associated  with  field  equation  solver. 

7.  Numerical  method  for  solving  the  field  equations 

Our  interest  in  this  work  was  to  simulate  the  nonlinear  wave  propagation  phenomena 
which  occur  in  the  high-speed  impact  of  munitions  on  targets.  Such  problems  involve 
highly  dynamic,  even  discontinuous  loading  of  the  system.  The  waves  generated  upon 
the  impact  loading  become  propagating  shock  waves.  Physically  realistic  weak 
solutions  of  the  governing  equations  are  therefore  sought.  There  are  a  host  of 
methods  designed  for  the  solution  of  hyperbolic  conservation  laws  in  the  presence  of 
shocks  developed  for  high-speed  aerodynamics  applications  (Leveque,  1990).  The 
challenge  in  this  report  is  to  adapt  these  methods  to  compute  wave  propagation 
phenomena  in  solids  where  the  constitutive  relationship  plays  a  role  in  the 
propagation  of  waves,  and  hence  deformation  of  the  media.  We  now  present  the  local 
Lax-Friedrichs  Essentially  Non-Oscillatory  (LLF-ENO)  scheme  used  for  solving  the 
conservation  laws.  In  order  to  apply  this  method  for  integration  of  the  equations,  it 
has  to  be  established  first  that  the  system  of  equations  under  consideration  is  indeed 
hyperbolic,  i.e.  that  the  eigenvalues  of  the  Jacobian  matrix  for  the  system  are  all  real. 
This  was  verified  to  be  the  case  for  the  range  of  physical  parameters  of  interest 
(material  properties,  velocities  etc.)  in  this  work  by  Vanden  (1998).  For  the 
homogeneous  form  of  the  ID  system  written  in  conservative  form,  the  eigenvalues 
are  found  to  be: 


A,  =  u 
-u 

Aj  =  M 


p 


These  eigenvalues  were  found  to  be  real  for  the  range  of  parameters  of  interest 
in  this  work. 


LLF-ENO  discretization  of  the  governing  equations: 


Consider  the  governing  equation  for  one-dimensional  transport: 


dt  dx 


Let 


dQ 

dt 


=  L(Q ) 


(21) 

(22) 
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Figure  5.  Illustration  of  grid  point  and  grid  face  definitions  for  discretization  of 
governing  equations.  H'+  and  FT  are  derivatives  of  the  interpolating  function  evaluated 
from  the  left  and  right  stencils  respectively. 


where 

L{Q)  =  -ZlILl  +  d{Q)  (23) 

*e-*w 

Fe  and  Fw  are  the  fluxes  at  the  east  and  west  faces  shown,  and  xe  and  xw  are  the 
locations  of  the  east  and  west  faces  respectively,  as  shown  in  Figure  5.  D  is  an 
appropriate  discrete  operator  for  the  source  terms.  In  the  current  work,  the  source 
terms  are  discretized  using  a  2nd-order  central  difference  scheme.  This  was  found  to 
be  robust  for  the  calculations  performed.  However,  it  may  be  necessary  in  future 
work  to  develop  a  more  sophisticated  differencing  procedure  for  the  source  terms  as 
well. 

The  three-step  third-order  in  time  Runge-Kutta  scheme  is  used  in  this  work  and  takes 
the  form: 


Qw  =  Qin)  +  AtL(Qin)) 

Qm  +  3Qin))  +  \AtL(Qil))  (24) 

4  4 

Qin+])  =^(2Qm  +QM)  +  ^AtL(Q{2)) 

The  spatial  order  of  accuracy  of  the  ENO  formulation  used  to  solve  Eq.  (1)  is 
determined  by  the  interpolation  practices  used  to  evaluate  the  fluxes  at  the  faces  e  and 
w,  i.e.  in  obtaining  Fe  and  Fw.  A  non-uniform  mesh  implementation  of  the  fluxes  in 
the  ENO  formulation  is  used  here  due  to  the  presence  of  immersed  boundaries,  as 
illustrated  in  Figure  6.  With  particular  reference  to  cell  j  ,  in  the  ID  case  the  interface 
can  lie  anywhere  between  Xj  and  Xj+i.  The  two  materials  are  treated  separately  in  the 
sharp  interface  formulation  and  the  flux  Fe  is  evaluated  at  a  location  (Xj  +  xim )  /  2 . 

Thus,  the  spacing  for  cell  j  is  different  for  the  e  and  w  faces.  The  flux  evaluations  for 
the  ENO  formulation  comes  from  derivatives  of  an  interpolating  function  H(x)  as 
follows: 
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Figure  6.  Grid  point  and  grid  face  definitions  for  evaluation  of  fluxes  in  the 
presence  of  an  immersed  boundary. 


F,  =£[//(*)  (25) 

The  derivatives  are  evaluated  from  divided  differences  and  the  flux  evaluation  is 
performed  as  follows: 


Fe  =  Fe+  +  Fe-  =H'+(xe)+H'-(xw)  (26) 

The  superscripts  (+)  and  (-)  indicate  the  positive  and  negative  direction  fluxes  at  the 
face  e  under  consideration  as  illustrated  in  Figure  5.  The  derivatives  H'  are  obtained 
as  explained  below.  Consider  the  interpolating  function  H(x).  In  terms  of  the  divided 
differences  this  function  can  be  written  as: 

H  (x)  =  H[xtl  ]  +  H[x„ ,  x]  ](x  -  x„ )  +  H[xa,  x, ,  x2  ](x  -  xa  )(x  -  x, )  +  0(h3 )  (27) 

his  interpolating  polynomial  can  be  carried  to  higher  orders.  We  will  restrict  attention 
here  to  developing  an  0(h2)  flux  approximation.  In  the  above  H[.,.]  symbolizes  the 
first  divided  difference  and  the  higher-order  divided  differences  are  obtained 
successively.  The  Essentially  Non  Oscillatory  (ENO)  schemes  (Shu  and  Osher, 
1988,1989),  are  derived  from  a  suitable  choice  of  the  stencil  locations  (x0,xi,X2,....) 
from  which  the  interpolating  function  is  constructed. 
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Figure  7.  Dlustration  of  stencils  used  for  obtaining  the  derivative  H’(+)  at 
face  Xj+i/2-  (a)  First-order  stencil,  (b)  Second-order  stencil  for  H’(2a)+. 

(c)  Second-order  stencil  for  H’(2b)+. 


For  example,  looking  at  Figure  7(a)  it  is  clear  that  there  is  only  one  stencil  possible 
for  the  first  divided  difference  for  control  point  location  j,  while  there  are  two 
candidate  stencils  (shown  in  Figures  7(b)  and  (c))  for  the  second  divided  difference, 
as  represented  by  the  forward  and  backward  differences  in  the  divided  difference 
table.  The  ENO  scheme  and  its  variants  derive  their  essentially  non-oscillatory 
property  from  the  choice  of  stencils  adopted.  The  original  ENO  scheme  (Harten  et  al., 
1987,  Shu  and  Osher,  1988,1989)  chooses  the  “smoothest”  stencil,  i.e.  the  lesser  of 
the  two  values  for  the  divided  differences  obtained  from  the  stencils  in  Figure  7(b) 
and  (c).  Weighted  ENO  schemes  (Jiang  and  Shu,  1996)  devise  appropriate  weights 
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for  each  candidate  divided  difference  and  then  evaluate  the  weighted  divided 
difference.  The  present  formulation  is  based  on  the  Convex  ENO  scheme  proposed  by 
Liu  and  Osher(1998)  and  chooses  the  divided  difference  value  “closest”  to  the 
previous  (here  the  first-)  order  flux  chosen.  It  is  this  adaptive  stencil  choice  procedure 
that  enables  the  Lax-Friedrichs  based  ENO  scheme  to  avoid  smearing  of  solutions 
away  from  shocks  while  maintaining  non-oscillatory  character.  Now  the  first  divided 
difference  is  obtained  as  follows: 


and 


H'+(x)  =  H< 


=  ^(f(q[Xj])+CC  Lq[xj } 


J+2 


H"(xe)  =  H- 


xe,x 


J+~2 


=  Lq[xj+l}> 

Z  J+  2 


Note  that  in  general: 

H +  [x0 ,  x,  }  ^  =  ^-(f(q(x))  +  a(x])q(x)) 
H~[x0,x^x=Xo  =^(f(q(x))-a(x0)q(x)) 

where 


-  (*0  +  *i) 

x  =  — - - — 


(28a) 


(28b) 


(29a) 

(29b) 


(30) 


These  of  course  apply  to  cells  away  from  the  immersed  boundary  such  as  cell  j+1  in 
Figure  6.  For  cells  that  are  adjacent  to  the  immersed  interface  such  as  cell  j  in  Figure 
6,  the  flux  evaluations  need  to  be  modified.  Here,  the  east  face  is  not  the  grid  cell  face 


but  is  located  at  ~ixj  +  xim )  where  Xjnt  is  the  location  of  the  interface.  Therefore,  for 


cell  y  : 


H'+(xe)  =  H[x  l,xiM]  =  Uf(q[xJ])+cc  iq[Xj}  (31a) 

j~2  2  J+2 

H'-(xe)  =  H[xe,xiM]  =  Uf(qim)-a  ,q-M)  (31b) 

2  7^ 


where  qin,  is  the  interfacial  value  of  the  convected  scalar  variable  q.  This  value  needs 
to  be  obtained  from  appropriate  boundary  conditions  applied  at  the  interface. 


The  first-order  flux  at  the  interface  is  then  obtained  using  Eqs  (26)  and  (28),  or 
Equations  (26)  and  (31)  if  cell  j  adjoins  the  interface.  In  order  to  determine  the  2nd- 
order  divided  difference,  the  following  steps  are  taken: 
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As  a  matter  of  notation  we  denote  the  first-order  flux  at  cell  face  located  at  x  ,  as  : 


f(j,j,+l)  =  H+[x  ,,x  ,  ]  =  ~-(f  (q(Xj))  +  a  xq[Xj}) 

1+~-,  2  J+x 


(32) 


where  the  notation  for  the  flux  at  the  face  in  terms  of  /(;1,;2,±1)  indicates  that  the 

flux  is  computed  for  the  face  of  cell  j2  using  values  at  control  point  jl.  The  ±1 
indicates  the  direction  of  flux  computed.  Therefore,  following  this  notation,  the  flux 
in  the  negative  direction  at  cell  face  located  at  x  ,  is  given  by 

j+l 

f(j  +  lj-l)  =  H~[x  ,,x  3]  =  ^-(f(q(x  ))-a  j q[Xj+] ])  (33) 

J+2  J+2  2  j+2 


The  candidate  second-order  derivatives  of  the  interpolating  function  H(x)  at  cell  face 
x  ,  are: 


H 


'(2  a)+ 


(*.  ,  )  =  tf  +  [x  3,X  ,  ]  + 


(x  3  +  x  ,  -2x  ,) 


J  ^  J  ~ 


j+~2 


/+2 


J-  j-  (x  |  —  X  3 ) 

J+2  i  2 


H  +  [x  3  ,  X  ,  ] 


.  3>* .  I 
J-2  J~2 


(2x  ,  -X  3  -X  ,) 

+  J2l  '~2  7~2 

(x  ,  -x  3) 

'+i  J~2 


H+[x  ,,x  ,] 

j~2  J+2 


(34) 


which  can  be  written  based  on  the  notation  in  Eqs.  (32)  and  (33)  as  : 


H'i2a)+  (x  +L )  =  f(j  - 1,  j,+ 1)  +  Aaa)+f(j  - 1,  j,+l)  +  Bi2a)+f(j,  j,+ 1)  (35) 

J  2 

where  A  and  B  are  grid  dependent  factors  determined  by  the  locations  of  the  stencil 
points  chosen.  Similarly,  the  other  candidate  second-order  fluxes  can  be  written  as: 


H'i2h)+(x  +, )  =  /  (j,  j,+ 1)  +  A(2h)+f(j,  j,+ 1)  +  Bm)+f(j  + 1,  ;,+l)  (36) 

J+2 

H'aa)~(x .+  , )  =  7 (j, 7,-1)  +  Ai2‘°-f(j , 7,-1)  +  B(2a)~f  (j  + 1, 7-I)  (37) 

J  2 

H'ah) ~  (x  +i )  =  f(j  + 1, 7,-1)  +  A(2b)~f  (7  + 1, 7,-1)  +  Bi2h)-f(j  +  2, 7,-1)  (38) 

1  2 
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Similar  expressions  for  the  3rd-  and  higher-order  derivatives  can  be  obtained. 
Therefore,  for  the  second-order  fluxes  at  the  face  j+1/2  there  are  two  candidates  each 
for  the  (+)  and  (-)  direction  contributions. 

In  the  presence  of  immersed  boundaries  the  discretization  in  the  cells  adjoining  the 
interface  only  will  need  to  be  modified  in  two  ways: 

1.  The  interface  boundary  conditions  will  appear  in  the  flux  contributions 
from  the  interface  side  as  in  Eq.  (31). 

2.  The  stencil  choices  possible  at  such  cells  will  be  limited  in  the  direction  in 
which  the  interface  lies.  For  example,  with  reference  to  Figure  6,  for  cell  j 
there  will  be  one  first-order  stencil  in  each  direction  as  for  the  interior 
cells.  However,  for  interfacial  cell  j+1,  there  can  be  only  one  choice  of 
second  order  flux  for  the  estimation  of  H'(  ) . 

These  considerations  are  no  different  in  fact  from  that  at  the  cell  immediately  in  the 
interior  of  the  domain  boundary.  Therefore  the  immersed  boundary  treatment  for 
evaluating  fluxes  is  no  different  from  that  for  domain  boundary  cells.  Apart  from 
these  considerations  the  fluxes  for  the  cell  j  adjoining  the  immersed  boundary  are 
constructed  using  Eq  (31)  except  that  the  values  of  qin) ,  i.e.  boundary  conditions  on 
the  immersed  interface  need  to  be  used  instead  of  the  grid  point  qj  values  in  the  Eq. 
(28).  It  is  not  readily  apparent  how  to  compute  the  boundary  values  for  all  the 
dependent  variables  in  the  particular  physical  problem  being  computed.  Some 
physically  based  boundary  conditions  can  be  imposed  based  on  the  physics  of  the 
impact  phenomena.  However,  some  of  the  physical  quantities  will  require  numerical 
boundary  conditions  as  in  other  systems  of  PDEs.  The  boundary  conditions  chosen 
and  the  rationale  for  the  choice  are  explained  in  the  following  sections. 


8.  Boundary  and  interface  conditions  in  the  ID  case 

If  the  rods  are  semi-infinite  the  boundary  conditions  at  x=0  and  x=L  ( L  being  the 
extent  of  the  domain  along  the  x-axis)  correspond  to  open  boundary  conditions  (i.e. 
zero-gradient  conditions).  However,  for  finite  size  rods  lying  within  the  domain  as 
shown  in  Figure  1 ,  two  types  of  boundary  conditions  need  to  be  considered. 


Type  1:  Material -material  (M-M)  boundary  condition: 

The  interface  between  the  two  rods  at  the  point  of  impact  is  one  separating 
two  materials.  When  the  two  materials  are  in  contact  continuity  dictates  that  the 
material  point  velocities  at  the  interface  are  equal,  the  total  stresses  are  equal  and  the 
temperatures  are  equal.  Therefore: 


=  u  =  umt 

(39) 

® x  ~~  ^*int 

(40) 

-T~  -T 

A  xint 

(41) 
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However,  there  are  six  dependent  variables  ( p,u,p,E,sx(ox),e )  in  the  governing 

equations  for  the  ID  system.  The  physically  available  boundary  conditions  and  the 
eos  yield  only  four  conditions  on  these  variables.  Therefore  one  needs  to  devise 
numerical  b.c.s  for  the  remaining  variables.  Several  alternative  formulations  of  b.c.s 
were  tried.  Some  of  these  yielded  unphysical  solutions.  The  viable  set  of  b.c.s 
determined  from  numerical  experiments  are  as  follows: 


P±  =  S(p7) 

(42) 

E+-=E(Ej) 

(43) 

c*=B(e,) 

(44) 

M*  =  I(Uj) 

(45) 

p±  =  eos(p±,(pu2)±,E±) 

(46) 

(Sx  +  p)±=I((sx  +  p)j) 

(47) 

Here  E  is  an  appropriate  extrapolation  operator.  The  choice  of  the  operator 
determines  the  order  of  accuracy  of  the  numerical  b.c.  at  the  interface.  We  perform 
linear  extrapolation  of  variables  from  the  suitably  chosen  grid  point  values  (e.g.  j+1 
and  j+2  on  the  (+)  and  j-1  and  j  on  the  (-)  side  of  the  interface  in  Figure  (6)  to  the 
interface.  Here  I  is  an  interpolation  operator.  Since  the  velocity  is  continuous  across 
the  interface  the  value  at  the  interface,  according  to  Eq.  (39)  can  be  obtained  by 
interpolating  from  the  grid  points  straddling  the  interface.  Again  the  order  of 
interpolating  operator  determines  the  accuracy  of  the  interface  b.c.  In  our  case  /  is  a 
distance-weighted  linear  interpolant. 


Type  2:  Material-void  (M-V)  boundary  condition: 

Before  the  rods  impact  there  is  a  void  separating  them.  Also,  if  the  rods  are  finite  the 
material  at  the  grid  points  to  the  right  of  the  right  rod  and  left  of  the  left  rod  are  void 
regions.  Therefore  the  ends  of  the  rods  are  material-void  interfaces.  Again,  the 
physics  dictates  the  following  b.c.s  on  the  material-void  interface: 

ax  =  <r+  =  o~  =  0  (48) 

along  with  the  equation  of  state  which  applies  in  the  material  at  the  interface.  This 
again  requires  imposition  of  numerical  boundary  conditions  at  such  interfaces  for  the 
remaining  variables.  The  appropriate  numerical  b.c.s  were  found  to  be: 


p~=E(Pj) 

(49) 

l 

II 

m 

(50) 

£~  =E(£j) 

(51) 
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(52) 

(53) 

(54) 


u~  =Z(Uj) 

p~=eos{p~,(pu2y,E~) 

(sx  +  pY  =  0 

Note  that  within  the  current  Eulerian-Lagrangian  framework,  the  boundary  conditions 
are  applied  and  the  discretization  is  performed  in  a  one-sided  manner,  i.e.,  the 
information  for  each  material  is  drawn  from  within  that  material.  Therefore,  at  the 
material-material  boundary  the  boundary  conditions  need  to  be  determined  on  both 
the  (+)  and  (-)  sides  of  the  interface  and  stored.  The  discretization  at  point  j  in  Figure 
(6)  therefore  uses  the  values  on  the  (-)  side  of  the  interface  while  that  at  point  j+1  uses 
the  values  on  the  (+)  side  of  the  interface.  Finally,  once  the  interface  velocity  is 
obtained  from  the  above  boundary  conditions,  the  interface  position  is  updated 
explicitly  as: 


x 


(H  +  l) 


=  X. 


(n) 


+  StU- 


(«) 


(55) 


9.  ID  Results 

The  ID  system  for  which  computations  were  performed  is  shown  in  Figure  1.  Since 
the  results  presented  below  are  for  the  case  of  uniaxial  strain,  the  two  materials  in 
impact  are  termed  plates.  The  plates  were  initially  placed  in  contact  at  the  midpoint  of 
the  domain  and  given  initial  velocities  which  caused  the  impact  and  the  interface 
separating  the  two  materials  then  traveled  over  the  fixed  grid.  The  position  of  both  the 
interface  separating  the  two  materials  as  well  as  the  ends  of  the  rods  (in  the  case  of 
finite  rods)  were  tracked  over  the  grid.  Thus  the  positions  of  the  plate  end  points  are 
tracked  explicitly. 

Results  of  three  cases  are  shown  in  the  following  for  the  ID  case.  These  involve  the 
traversal  of  plate  end  points  of  the  material-material  and  material-void  interfaces  on 
the  fixed  mesh.  All  computations  are  performed  using  the  material  properties  for 
Copper.  The  values  used  were  a  shear  modulus  of  G=4.5xl010  Pa,  density 
po=8.93xl03  kg/m3,  and  an  idealized  stress-strain  relationship  giving  the  equivalent 
stress: 

at  =  Jfyij  =  #  (£)  =  H0  +  H'e  (56) 

where  H0=1.2OxlO9  Pa  and  a  constant  H'=  0.12xl09Pa.  Also  in  the  Mie-Gruneisen 
equation  of  state  for  copper,  To=2,  co=3.94xl03  m/s  and  s=1.49. 

The  mesh  used  in  each  case  consists  of  100  points  in  the  domain.  There  are  four 
interfacial  points  tracked  on  the  mesh  corresponding  to  the  ends  of  the  two  bars.  The 
computations  are  carried  to  significant  times  after  impact  so  that  the  interfaces  move 
over  several  mesh  cells. 
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In  Figure  8,  we  show  the  case  of  impact  of  two  plates  of  semi-infinite  extent.  The 
plates  are  initially  positioned  such  that  impact  takes  place  at  x  =  0.6x1  O'5,  the  center 
of  the  domain.  The  left  plate  has  a  velocity  before  impact  of  500  m/s.  The  interface 
points  (plate  end  points)  at  the  center  in  the  domain  are  supplied  with  the  material- 
material  boundary  conditions  on  the  (+)  and  (-)  sides.  At  the  left  and  right  end  of  the 
domain,  zero-gradient  conditions  are  applied  for  all  quantities  since  the  plate  is 
assumed  to  be  semi-infinite.  In  each  figure  we  show  the  profile  along  x  of  the  flow 
variables  at  equal  intervals  of  time  (At=2  microseconds)  after  impact.  In  figure  8(a) 
the  velocity  profiles  are  shown.  In  Figure  8(c)  we  show  the  phases  at  each  point  in  the 
grid,  again  at  equal  time  intervals  of  2  microseconds.  It  can  be  seen  that  the  interface 
travels  to  the  right  due  to  the  impact.  In  Figure  8(b)  the  profiles  of  the  equivalent 
stress  show  the  large  gradients  developed  at  the  interface  following  impact.  The 
resulting  wave  propagation  phenomena  lead  to  the  propagation  of  stress  waves  and 
corresponding  strains  through  the  plates.  Symmetry  about  the  point  of  impact  is 
maintained  during  the  computation  as  it  should  be.  At  the  material-material  interface 
point  the  material-material  boundary  conditions  discussed  above  have  been  used. 
Clearly,  in  the  presence  of  the  moving  material-material  boundary,  the  solution  of  the 
governing  equations  leads  to  a  stable  physical  solution. 

In  Figure  9  we  show  the  solutions  for  the  case  where  there  are  two  finite  plates  in  the 
domain.  The  left  plate  initially  has  a  velocity  to  the  right  of  500m/s.  The  right  plate 
is  initially  stationary.  The  finite  plates  initially  occupy  the  regions 
lxlO-5  <  x,  <  3xlO-5  for  the  left  plate  and  3xl0~5  <  jcr  <  5jc10-5  for  the  right  plate. 
We  show  the  profiles  for  velocity  (Figure  9(a)),  equivalent  stress  (Figure  9(b)),  and 
phases  (Figure  9(c))  in  the  domain  at  equal  intervals  of  time  of  At=2  microseconds. 
In  this  case,  while  the  waves  propagate  in  the  two  plates,  the  plates  are  physically 
moving  and  deforming  over  the  mesh.  The  gradual  equalization  of  the  velocities  of 
the  material  points  in  the  two  plates  is  seen  in  Figure  9(a).  As  the  computation 
proceeds  the  maximum  material  point  velocity  oscillates  between  the  two  plates.  As 
can  be  seen  in  the  profile  7  in  Figure  9(a),  the  particle  velocity  in  each  plate  is  tending 
to  the  mean  value  of  250m/s.  Immediately  after  impact  we  see  very  large  stresses 
developing  at  the  point  of  impact.  As  can  be  seen  in  Figure  9(c)  the  two  plates  are 
moving  to  the  right  and  the  four  interfaces,  namely  the  material-material  and 
material-void  interfaces  are  tracked  through  the  mesh. 
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(a) 


(c) 


(b) 


Figure  8.  ID  calculation  of  the  impact  of 
two  slender  infinitely  long  plates  on  the 
x-axis.  Profiles  of  flow  quantities  at 
equal  intervals  of  time  (2  microseconds) 
are  shown: 

(a)  velocity  profiles  (b)  Stress  (c)  phase 
(phase=250  corresponds  to  left  plate, 
phase=500  corresponds  to  right  plate). 
Numbers  on  the  profiles  indicate  the 
sequence  in  time. 


In  Figure  10,  we  compute  the  case  of  two  finite  plates  impacting  with  opposite 
velocities.  Initially  the  plates  occupied  the  regions,  lx  10'5  <  x,  <  6xl0”5,  and 

6xl0~5  <  xr  <  l.lxlO-5.  The  left  plate  was  given  a  velocity  of  250m/s  and  the  right 
plate  of  -250  m/s.  Note  in  this  case  the  final  velocities  by  symmetry  should  tend  to 
cancel,  and  both  plates  should  be  stationary  on  the  x-axis,  i.e.  one  expects  that  the 
wave  propagation  and  reflection  from  the  material-void  interfaces  at  the  far  ends  of 
the  plate  should  take  place  symmetrically.  Clearly,  this  symmetry  is  well  maintained 
by  the  computations.  Initially,  as  seen  in  Figure  10(b)  large  gradients  in  stress  are 
initially  developed  in  the  region  close  to  impact.  The  longitudinal  oscillation  of  the 
plates  about  the  mean  position  can  be  seen  from  the  profiles  in  Figure  10. 
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(a) 


(c) 


(b) 


Figure  9.  ID  calculation  of  impact  of  two  slender 
finite  length  plates  on  the  x-axis.  Profiles  of  flow 
quantities  at  equal  intervals  of  time  (2 
microseconds)  are  shown: 

(a)  velocity  profiles  (b)  Stress  (c)  phase 
(phase=250  corresponds  to  left  plate,  phase=500 
corresponds  to  right  plate).  Numbers  on  the 
profiles  indicate  the  sequence  in  time. 


10.  Boundary  and  initial  conditions  in  the  2D  case 


A  schematic  for  the  interface  treatment  in  2D  computations  is  shown  in  Figure  1 1 . 
The  impactor  and  target  are  initially  placed  some  distance  apart  and  impact  is 
initiated  by  giving  a  uniform  velocity  to  the  impactor  while  the  target  remains 
stationary.  Therefore,  initially  only  particle  velocities  and  other  material  properties 
such  as  density  have  non-zero  values,  while  stresses  are  set  to  zero.  In  the  2D  case  the 
boundary  conditions  are  based  on  the  physically  imposed  ones  and  numerical  b.c.s 
obtained  based  on  the  experience  gained  from  the  ID  situation  and  tested  by 
numerical  experiments.  In  contrast  with  the  ID  case  however,  in  2D  some  regions  of 
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(a) 


(b) 


(c) 


Figure  10.  ID  calculation  of  impact  of 
two  slender  finite  length  plates  on  the  x- 
axis.  Profiles  of  flow  quantities  at  equal 
intervals  of  time  (2  microseconds)  are 
shown: 

(a)  velocity  profiles  (b)  Stress  (c)  phase 
(phase=250  corresponds  to  left  plate, 
phase=500  corresponds  to  right  plate). 
Numbers  on  the  profiles  indicate  the 
sequence  in  time. 


the  interface  can  be  in  material-material  contact  and  others  can  be  exposed  to  a  void. 
This  is  illustrated  in  Figure  11.  Here  the  region  of  contact  between  two  curved 
boundaries  in  collision  is  shown.  The  individual  marker  points  which  define  the  two 
interfaces  are  also  shown  in  the  Figure  11.  Interfacial  boundary  conditions  are  applied 
at  these  marker  locations  on  the  interface.  Therefore,  as  can  be  seen  in  the  figure,  an 
interfacial  marker  in  the  material-material  region  can  have  an  immediate  neighbour  in 
the  material-void  region  of  the  interface.  Therefore,  at  each  instant  of  the  interface 
deformation  the  interface  markers  have  to  be  classified  as  material-material  or 
material-void  markers  and  the  appropriate  b.c.s  obtained.  Thus  two  immediately 
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Figure  11.  Illustration  of  marker  location  and  the  type  of  boundary  condition 
applied  at  the  different  interface  regions. 


adjacent  interfacial  marker  neighbours  can  have  entirely  different  b.c.s  imposed  on 
them. 


Type  1:  Material -material  interface 

Boundary  conditions  in  the  material-material  case  are  developed  based  on  the 
physically  required  conditions  of  equal  material  point  velocities  at  the  interface  for 
the  two  materials,  and  the  continuity  of  stress  and  temperature.  These  conditions  can 
be  stated  as: 


u+=U-=  wint 

(57) 

II 

1 

II 

3 

(58) 

^  xx  ^  xx  (^jor)int 

(59) 

<7,v+  =  Cy>:  =  (<T,v)int 

(60) 

<*xy  —  ^ xx  ~  (^jry)int 

(61) 

T+=T-=Tim 

(62) 

Along  with  the  eos  for  the  pressure  these  constitute  7  physically  imposed  conditions 
at  the  interface  for  the  nine  independent  variables  ( p,u,v, p,E,sxx,s  ,s  ,£).  This 

necessitates  numerical  boundary  conditions  to  be  developed  for  the  remaining 
dependent  variables.  We  obtain  these  based  on  extensions  of  the  ID  boundary 
conditions  and  verify  that  they  provide  numerically  stable,  convergent  solutions 
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Figure  12.  Illustration  of  marker  and  grid  points  involved  in  applying  the  material- 
material  boundary  condition. 


without  unphysical  behaviors  such  as  overshoots  or  oscillations.  The  numerical 
implementation  of  the  boundary  conditions  can  be  explained  with  reference  to 
Figures  (12)  and  (13). 

For  the  material-material  contact  situation  shown  in  Figure  12,  the  boundary 
conditions  are  imposed  on  the  interfacial  markers  shown  on  the  interface  by  the  filled 
circles.  The  viable  boundary  conditions  are  determined  to  be: 


P±=S(P,7) 

(63) 

E±=E(EiJ) 

(64) 

£+-=m,j) 

(65) 

p±  =eos(p±,(pu2)±,(pv2)±,E±) 

(66) 

«*  =  l(ui}) 

(67) 

v±=/(v,P 

(68) 

{sxx  +  pf  =  l{(sxx  +  p).j) 

(69) 

(syy  +  p)±  =  I((Syy  +  p)JJ) 

(70) 

(^v)±=/((^v)..) 

(71) 

In  the  above  the  operators  E  and  I  are  the  extrapolation  and  interpolation  operators 
respectively.  The  subscripts  (i,j)  indicate  values  on  the  grid  while  superscripts  ± 
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Void 


Figure  13.  Illustration  of  marker  and  grid  points  involved  in  applying  the  material-void 
boundary  condition. 


indicate  values  on  the  (+)  and  (-)  sides  of  the  interface.  In  particular  for  the 
interfacial  marker  indicated  by  the  ellipse  in  Figure  12,  a  bilinear  interpolation 
operator  is  used  to  estimate  the  variables  which  are  continuous  across  the  interface. 
The  points  1,2,3  and  4  straddle  the  interface  marker  in  question  and  bilinear 
interpolation  is  performed  based  on  the  location  of  the  interface  marker  and  the 
location  of  these  grid  points.  Note  that  on  a  Cartesian  mesh,  having  obtained  the 
information  regarding  the  interface  and  its  relation  to  the  mesh  as  discussed  in 
Section  6,  it  is  a  simple  matter  to  choose  the  points  in  the  bilinear  interpolant.  The 
extrapolation  operator  used  for  obtaining  the  values  of  variables  not  governed  by 
continuity  conditions  is  somewhat  less  straightforward  in  2D.  For  example,  with 
particular  reference  to  the  interfacial  marker  in  focus  in  Figure  13  (in  the  ellipse), 
extrapolation  has  to  be  performed  from  values  at  grid  points  within  the  interior  of  the 
material.  In  order  to  do  so,  a  normal  is  extended  from  the  interfacial  marker  in 
question  into  the  material.  The  values  at  two  nodes  on  the  normal  placed  a  distance  h 
(the  grid  spacing)  apart  are  obtained  by  bilinear  interpolation  from  grid  points 
straddling  each  normal  node  point.  For  normal  node  point  1  for  example,  the  values 
are  interpolated  from  grid  points  connected  by  dashed  lines.  The  value  at  the  interface 
is  then  obtained  by  extrapolating  from  the  values  at  the  two  nodes  along  the  extended 
normal  to  the  interface. 

Type  2:  Material-void  interface 

For  a  material-void  interface,  the  physically  imposed  conditions  on  the  interface  are 
that  the  traction  be  zero.  In  the  present  case  we  let  the  stresses  vanish  on  the  free 
surface  as: 
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(72) 

(73) 

(74) 


Along  with  the  eos  for  pressure,  which  applies  at  the  interface  on  the  material  side, 
the  imposed  b.c.s  account  for  only  four  of  the  nine  independent  variables.  Therefore 
numerical  b.c.s  are  required  at  the  interface  and  these  are  devised  based  on  extensions 
from  the  ID  case  and  by  experimentation  such  that  the  b.c.s  imposed  do  not  lead  to 
development  of  unphysical  features  in  the  flowfield.  The  b.c.s  imposed  on  the 
material  side  of  the  interface  are  as  follows: 


P"  =  S(p„,.) 

(75) 

iT=E(£,.) 

(76) 

e-=E(e,j) 

(77) 

u  =  5  (Uij) 

(78) 

v-  =  S(v. ,) 

(79) 

P "  =eos(p~,(pu1y,(pv2y,E~ ) 

(80) 

+ 

1 

II 

o 

(81) 

3* 

+ 

1 

II 

o 

(82) 

(O"=0 

(83) 

Note  that  no  b.c.s  are  required  on  the  void  side  of  the  interface. 


11.  Results  of  2D  computations 

2D  computations  were  performed  in  a  square  domain  of  size  lm  x  lm  as  illustrated 
in  Figure  2.  As  shown  there  the  objects  were  placed  some  distance  apart  on  the  mesh 
and  impact  was  initiated  by  prescribing  a  velocity  to  one  or  both  interfaces. 
Initially  there  is  a  region  of  void  between  the  two  interfaces.  When  the  interface 
markers  lie  in  a  material-void  region  of  the  interface  we  apply  the  boundary 
conditions  of  type  2  (M-V)  on  such  markers.  On  markers  in  the  material-material 
regions  of  the  interface,  b.c.s  of  type  1  (M-M)  are  applied.  Note  that  the 
discretization  procedure  in  2D  requires  the  value  of  qint  at  interface  location  not 
coincident  with  the  interfacial  markers.  To  obtain  the  interfacial  value  of  flow 
quantities  at  such  points,  we  interpolate  along  the  interface  using  a  linear  function  fit 
between  successive  markers  straddling  the  point  as  shown  in  Figure  4(c).  Higher- 
order  interpolation  along  the  interface  is  of  course  possible  but  was  found  to  lead  to 
erroneous  behaviour  due  to  the  large  gradients  along  the  interface  at  regions  where 
there  is  an  abrupt  switch  from  a  M-M  to  M-V  type  boundary  condition.  This  implies 
that  at  points  where  the  interfaces  go  from  being  in  the  material-material  contact 
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situation  to  material-void  situation,  such  as  point  A  in  figure  11,  we  introduce  a 
lower-order  interpolation  practice  to  avoid  overshoots. 

We  have  tried  several  challenging  cases  for  the  evolution  of  the  interfaces  after 
impact.  In  the  first  case,  presented  in  Figure  14,  we  show  the  impact  of  a  copper 
cylinder  with  a  copper  planar  surface.  Both  surfaces  are  copper  and  the  material 
properties  in  the  model  correspond  to  that  metal.  In  the  figure,  we  show  on  the  left  the 
contours  of  particle  velocity  magnitude  in  the  impactor  and  the  target  along  with  the 
velocity  vectors  in  the  flow  domain.  On  the  right  we  show  contours  of  equivalent 
stress.  Also  shown  in  each  of  the  figures  is  the  shape  of  the  boundaries  of  the  two 
materials.  As  can  be  seen  in  these  figures  there  is  an  abrupt  transition  in  the  corners 
from  a  mterial-material  interface  to  a  material-void  interface  for  each  material. 
Appropriate  boundary  conditions  as  discussed  in  Section  10  are  applied  in  these 
regions.  Zero-gradient  conditions  are  applied  at  the  sides  of  the  domain  assuming 
that  the  target  has  infinite  extent  in  all  except  the  +y  direction.  Figures  14(a),  (b)  and 
(c)  correspond  to  time  instants  2.5  ps,  50|is  and  100  ps  after  impact  respectively.  The 
progression  of  the  elastic-plastic  waves  and  the  formation  of  large  gradients  in  the 
velocity  as  well  stress  fields  is  evident  from  the  figure.  At  the  rim  of  the  impactor, 
the  interfaces  are  constantly  in  collision  since  the  material-void  interfaces  are  being 
pushed  against  each  other  to  form  material-material  interfaces.  Thefore  the  rim  of  the 
impact  region  registers  large  stress  and  correspondingly,  strain  values.  Stress  waves 
are  propagated  into  the  materials  from  this  point.  In  Figure  14(c)  it  can  be  seen  that 
the  velocity  field  is  such  as  to  continuously  push  the  impactor  into  the  target  leading 
to  the  production  of  an  upswell  in  the  target  material  around  the  rim.  This  is  also 
indicated  clearly  by  the  velocity  vectors  shown.  Regions  of  compression  and  tension 
are  seen  from  the  contours  of  stress.  The  computational  time  for  a  100x100  mesh 
calculation  to  the  stage  shown  in  Figure  14(c)  is  about  one  hour  on  a  DEC- Alpha 
5000  computer.  No  attempt  has  been  made  to  optimize  the  code  and  it  should  be 
possible  to  decrease  computational  times  significantly. 

The  next  case  corresponds  to  the  impact  of  two  cylindrical  objects.  Initially  the  lower 
object  is  stationary  and  the  upper  object  is  placed  some  distance  away.  The  upper 
object  then  is  projected  with  a  velocity  of  -2000  m/s  and  impact  with  the  lower 
object  results  in  the  situation  shown  in  the  figure.  Figures  15(a)  ,(b),(c)  and  (d) 
correspond  to  2.5  ps,  50  ps,  100  ps  and  150ps  after  impact  respectively.  The  lower 
object  acquires  velocity  after  impact  and  is  allowed  to  move  freely  downward 
through  the  mesh  as  the  velocity  field  develops.  The  velocity  vectors  show  that  the 
impact  causes  a  flattening  of  the  interfaces  and  that  the  lower  target  deforms 
symmetrically  with  respect  to  the  impactor.  Again  the  rim  of  the  impact  region 
shows  large  stress  concentrations  due  to  the  continuous  compression  of  the  interfaces 
at  that  point  by  incorporation  of  the  material-void  interfaces. 
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(a)  t=2.5  microseconds 


(b)  t=50  microseconds 


wmm 

StL 

m 

I 


(c)  t=  100  microseconds 


Figure  14.  Impact  of  a  cylinder  with  a  planar  surface.  The  cylinder  impacts  the  target  with  a 
velocity  of  2000m/s  directed  downward.  The  figures  on  left  show  velocity  contours  and 
vectors  along  with  the  interface  shapes.  The  times  after  impact  are  indicated  alongside  the 
figures.  The  figures  on  right  show  stress  contours. 
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Figure  15 


(d)  t=150  microseconds 


Figure  15  (continued).  Impact  of  two  cylinders  in  the  plane.  The  upper  cylinder  impacts 
the  target  with  a  velocity  of  2000m/s  directed  downward.  The  figures  on  left  show 
velocity  contours  and  vectors  along  with  the  interface  shapes.  The  times  after  impact  are 
indicated  alongside  the  figures.  The  figures  on  right  show  stress  contours. 


In  Figure  16  we  show  the  effect  of  grid  refinement  and  time  step  refinement  on  the 
computation  of  this  problem.  In  that  figure  we  have  plotted  the  interfaces  computed 
at50ps,  75ps,  100  ps  and  15 Ops  respectively  for  the  following: 

1 .  50x50  mesh,  CFL  =0.5 

2.  100x100  mesh,  CFL  =0.5 

3.  100x100  mesh  ,  CFL  =0.1 

It  can  be  seen  from  the  figure  that  the  solutions  for  these  three  cases  are  close 
throughout  the  calculation.  In  fact  cases  2  and  3  above  are  almost  indistinuishable 
due  to  the  high-order  time  integration  used  (  3rd-order  ).  However,  for  t=100  ps  and 
150  ps  one  sees  differences  in  the  solutions  for  the  two  grid  resolutions.  The  solution 
for  the  coarser  grid  lags  that  of  the  finer  grid.  Also  the  largest  differences  appear  to 
be  at  the  rims  of  the  impact  area  where  the  largest  gradients  are  encountered  due  to 
the  singular  behaviour  there. 

As  a  final  demonstration  we  simulate  the  impact  of  a  projectile-like  slender  body  with 
a  planar  target.  The  results  are  shown  in  Figure  17.  The  features  of  the  impact 
process  in  this  case  are  very  similar  to  those  observed  before.  One  again  sees  large 
stress  concentrations  in  the  rim  of  the  impact  zone  and  the  upwelling  of  material  of 
the  target  in  the  region  around  the  impact  zone.  This  is  also  seen  from  the  velocity 
vector  plots.  The  direction  of  particle  velocities  reverse  in  the  impacted  material 
away  from  the  impacted  region  leading  to  the  upwelling  of  material.  This  is  despite 
the  specification  of  zero  gradient  conditions  at  all  boundaries  in  the  target  which 
were  specified  to  approximate  the  infinite  extent  of  the  target.  The  presence  of  the 
free  surface  (material-void  interface)  surrounding  the  impacted  region  allows  for  the 
development  of  a  velocity  at  the  free  surface  which  perits  the  pushing  and  lifting  of 
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Figure  16  Grid  refinement  study  for  cylinders  impacting  in  the  plane.  The 
interface  shapes  are  superposed  for  the  following  cases:  50x50  mesh,  CFL=0.5; 
100x100  mesh,  CFL=0.5;  100x100  mesh,  CFL=0.1.  Differences  can  only  be 
discerned  between  the  50x50  and  100x100  mesh  solutions.  Interfaces  for  the 
two  different  time  step  sizes  lie  on  each  other. 


material  in  that  region.  Note  the  shape  of  the  crater  formed  due  to  the  impact.  Since 
the  target  is  unconfined  the  resulting  crater  appears  to  be  assuming  a  nearly  circular 
shape. 

11.  Summary  and  future  work 

We  have  described  the  development  of  a  numerical  technique  based  on  a  fixed-grid 
sharp  interface  tracking  approach  for  the  simulation  of  multi-material  impact.  The 
physics  of  the  problem  is  such  that  nonlinear  elastic-plastic  wave  propagation 
phenomena  occur  in  the  materials  leading  to  the  formation  of  shocks.  We  track 
interfaces  explicitly  as  points  in  ID  and  curves  in  2D.  In  its  interaction  with  the 
flowfield,  although  we  are  computing  on  a  fixed  mesh,  the  interface  is  treated  sharply 
and  the  discontinuities  at  the  interface  are  not  smeared.  We  have  demonstrated  that 
the  current  method  has  the  following  capabilities: 
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(d)  t=  1 50  microseconds 


Figure  17.  (Continued)  Impact  of  a  slender  elliptic  cylinder  with  a  planar  surface.  The 
cylinder  impacts  the  target  with  a  velocity  of  2000m/s  directed  downward.  The  figures  on  left 
show  velocity  contours  and  vectors  along  with  the  interface  shapes.  The  times  after  impact 
are  indicated  alongside  the  figures.  The  figures  on  right  show  stress  contours. 


1 .  The  interface  can  be  tracked  through  large  distortions. 

2.  Accurate  shock-capturing  schemes  can  be  implemented  for  Cartesian  grids 
and  extended  in  a  straightforward  manner  to  incorporate  the  presence  of 
the  moving  interfaces. 

3.  Boundary  conditions  are  developed  for  the  ID  uniaxial  strain  case  and  2D 
plane  strain  case  and  these  are  applied  at  the  exact  locations  of  the 
boundaries. 

4.  Different  regions  of  the  boundaries  can  have  different  boundary 
conditions,  i.e.  the  material-material  and  material-void  boundary 
conditions.  These  are  applied  at  the  interface  points  identified  to  lie  in 
regions  where  the  interfaces  are  in  contact  and  where  the  interface  is 
exposed  to  a  void,  respectively.  These  boundary  conditions  are  a 
combination  of  physical  and  numerical  boundary  conditions.  The 
suitability  of  the  set  of  b.c.s  is  determined  based  on  numerical 
experimentation.  The  singularity  resulting  from  an  abrupt  transition  from  a 
material-material  to  material-void  b.c.  at  the  interfaces  is  handled  well. 

5.  Computations  of  the  deformation  process  are  carried  to  large  distortions 
while  the  interfaces  travel  through  the  mesh  in  a  stable  and  robust  manner. 
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Future  work  will  be  directed  toward  establishing  the  validity  of  this  method  by 
solving  problems  that  can  be  benchmarked,  such  as  the  Taylor  bar  impact  problem  in 
2D.  We  will  also  work  on  augmenting  the  current  capabilities  by  introducing  grid 
refinement  adaptively  in  the  vicinity  of  interfaces  and  regions  of  large  gradient.  This 
will  alleviate  computational  times  and  also  enhance  accuracy  in  desired  regions  thus 
enabling  large-scale  computations. 
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